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ABSTRACT 

We have adapted Coupled Escape Probability, a new exact method of solving radiative transfer 
problems, for use in asymmetrical spherical situations. Our model is intended specifically for use in 
modeling optically thick cometary comae, although not limited to such use. This method enables 
the accurate modeling of comets' spectra even in the potentially optically thick regions nearest the 
nucleus, such as those seen in Deep Impact observations of 9P/Tempel 1 and EPOXI observations of 
103P/Hartley 2. 

Subject headings: Comets, Methods: numerical, Opacity, Radiative Transfer 



1. INTRODUCTION 

Comets are understood to be frozen remnants from the 
formation of our solar system. As such, their chemical 
composition is of great significance to understanding the 
origin of the planets and the distribution of important 
molecules, including water, throughout the solar system. 
This was and is a major goal of the Deep Impact and 
EPOXI Missions, among others, as well as ground-based 
observations of comets. 

Recent observations, in particular those of the Deep 
Impact and EPOXI Missions (see e.g. Feaga, et al. 2007 
or A'Hearn, et al. 2011), have provided better spectra 
of a cometary coma than were, in general, previously 
available. These observations include spectra with high 
spatial resolution very near to the nucleus. No previous 
observations had as much spatially well-resolved spec- 
tral data, and thus there was little observationally driven 
need to pay special or close attention to the densest part 
of the coma. Ground-based observations could only see 
optically thick regions of comae for the brightest and/or 
most active of comets, (e.g. Hale Bopp; see DiSanti et 
al. 2001.) 

Therefore many earlier studies that modeled spectra of 
comae, in keeping with the available observations of the 
time, did not attempt to calculate optical depth effects 
on spectra. Optically thin comae were assumed, since the 
field of view in those observations being modeled would 
be dominated by the majority of the coma far from the 
nucleus, which is optically thin (see e.g. Chin & Weaver 
1984; Crovisier & Le Bourlot 1983). However, with the 
proliferation of space missions to comets, as well as much 
better instruments for ground-based observations (see, 
e.g. DiSanti et al. 1999), this is no longer a truly tenable 
approach. 

Our goal is to better understand the abundances, dis- 
tributions and creation mechanisms of various gases ob- 
served in comae, in particular of comet 9P/Tempel 1, the 
target of the Deep Impact Mission, and 103P/Hartley 2, 

agcrsch@astro.umd.edu 



the subject of the EPOXI mission. 

In order to do so, we have built a computer model of 
the spectrum of the comet's coma which includes the dif- 
ficult and often ignored problem of accurately including 
radiative transfer to account for the potentially optically 
thick coma (or regions of the coma) near the nucleus. 

This model will facilitate analyzing the actual spectral 
data from the Deep Impact and EPOXI missions to bet- 
ter determine abundances of key species, including CO, 
C0 2 , and H 2 0, as well as remote sensing data on active 
comets. 

1.1. The Simple Coma Model 

We begin our modeling of IR ro- vibrational spectra of 
a coma by initially following the method used by Chin 
& Weaver 1984, Crovisier 1987, and others with some 
minor improvements. As in those models, we assume 
a constant expansion velocity, thus linearly relating any 
radial distance to a specific time since a "parcel" of gas 
was released from the surface of the nucleus. Therefore 
we numerically integrate over time the linear differential 
equations defined by the Einstein coefficients and colli- 
sional rate coefficients to get fractional molecular enegry 
level populations for each distance from the nucleus, from 
which we calculate emission spectra. 

^ = 5> 4fe n, + J Vik {B ikni - B kl n k )) + (1) 
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Here i, j, k and / indicate energy levels of a molecule with 
the n's with those indices being the corresponding level 
populations. The Einstein coefficients between levels x 
and y are A xy and B xyi and C is a similar collisional 
coefficient. We use J„ for the mean intensity of radi- 
ation at the frequency corresponding to the transition 
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between x and y, E xy is the energy difference between 
levels x and y, and fee and T have their usual meanings. 
The summations are over all levels i, 1 or j which have 
a transition into or out of level k. For collisional coeffi- 
cients C — n H ^ Q av where n H ^ Q is the number density of 
H 2 (assumed to be the dominant collisional partner) a 
is the collisional cross section for a given transition and v 
is the mean (thermal) molecular speed. This allows us to 
include a time variable production rate. We use such a 
coma integration as the initial basis for our coma model 
before including potential optical depth effects. 

Our primary improvement is the inclusion of radia- 
tive transfer calculations using our spherical adaptation 
of the Coupled Escape Probability method (hereafter, 
"CEP"; see Elitzur & Asensio Ramos 2006, hereafter, 
"CEP06") to more correctly model optically thick (or 
potentially thick) regions of cometary comae. This is de- 
scribed in detail below, and is the main part of this paper. 
We use the coma integration results to provide the "ini- 
tial guess" values for populations used in the subsequent 
radiative transfer calculations using CEP. 

For the purposes of the initial coma model, we treat 
the comet as spherically symmetric, and as having a uni- 
form and constant gas production rate over its entire 
surface. The outward speed of the gas is also assumed 
to be constant, as per Chin & Weaver 1984. While this 
is not strictly physically accurate (see e.g. Combi 1996) 
the variation over the majority of the coma is relatively 
small. We use a temperature profile that varies with ra- 
dial distance from the nucleus, having closely followed 
Combi's 1989 model (see Fig. [I}) These approxima- 
tions make integration over time equivalent to calculat- 
ing these values over increasing distances from the comet 
nucleus for a "shell" of gas expanding outwards from the 
nucleus. 
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Fig. 1. — Temperature profile of cometary atmosphere based on 
Combi 1989 for surface temperature of 200 K and radial gas velocity 
of 0.8 km s" 1 . 

We ignore the photodestruction of CO in our coma 
model, due to its long lifetime (see Crovisier 1994.) The 
lifetime is > 10 7 s, and we are integrating out to 10 5 km 
with an expansion velocity of 0.8 km s _1 . (Note that 
others, e.g. Morgenthaler, et al. 2011, find a shorter 
lifetime, but still >4x 10 5 s, which is still large enough 
that it can be neglected in our modeling out to 10 5 km.) 

The model can include coma morphology features as 
well, each modeled with its own coma integration using 
separate conditions. Such features, as seen in the Deep 
Impact and EPOXI encounters, are a main motivation 
for creating this model to better understand possible op- 
tical depth effects in the near nucleus regions of the coma. 
(See e.g. Feaga, et al. 2007 or A'Hearn, et al. 2011) 



After describing our method in Section 2, we present 
(in Section 3) some results demonstrating its use in bet- 
ter understanding possible optically thick spectra for the 
carbon monoxide 1-0 (X 1 !]"*") band in spherically sym- 
metric comae. These may be useful for ground-based 
observers (or space telescope observations of comets) to 
better fathom the depths of cometary spectra. Forth- 
coming model results for CO2, H 2 and near- nucleus 
morphological features will follow (in other papers). 

2. OUR METHOD: CEP ADAPTED FOR AN ASYMMETRIC 
SPHERICAL CASE 

This section describes our adaption of the Coupled Es- 
cape Probability radiative transfer technique to spherical 
cases in which the plane parallel approximation is not ap- 
propriate, including most cometary problems. (Note that 
Yun, et al. 2009, previously adapted CEP for purely sym- 
metric spherical cases.) We have developed this model 
specifically for use in studying cometary comae, but it 
could be applied to many other astrophysical phenomena 
as well (e.g. planetary atmospheres, molecular clouds, 
etc.). 

We derive expressions analogous to the relevant CEP06 
equations, adapted from the plane parallel situation to 
spherically symmetric cases. Then we include asymme- 
tries of two different types: radiation from an outside 
source (here, the Sun) and non-uniform conditions due 
to morphology. 

In brief, the CEP method (see CEP06) divides up a 
plane-parallel "slab" into "zones" (each of which has uni- 
form properties) and calculates the net radiative bracket 
("p") that is multiplied by the Einstein A coefficient in 
the equations of statistical equilibrium (see Equation [IJ 
for each radiative transition for each zone. (Note that 
the "p" term effectively combines the "B" terms into the 
"A" term. See CEP06 for more detail.) The innovation 
of CEP is that the net radiative brackets for each zone 
can accurately represent the contributions of all zones' 
emission and absorption to/from other zones. (As op- 
posed to "plain" Escape Probability where a similar fac- 
tor added to the statistical equilibrium equations is only 
a local approximation of a photon's likelihood of escaping 
the entire slab. See, Bockelee-Morvan 1987, Zakharov, 
et al. 2007 and Litvak & Kuiper 1982.) The statistical 
equilibrium equations for all zones, with the inclusion of 
the net radiative bracket, form a single non-linear ma- 
trix. This matrix can be solved using an algorithm for 
non-linear matrix solving such as Newton's Method. We 
use functions from Numerical Recipes in C (see Press, 
et al. 1992) for the Newton based matrix solver, as well 
as other calculations such as numerical integration. Fur- 
ther computational details can be found in the Appendix. 
This solution yields the fractional populations of molec- 
ular energy levels for each zone. From these, the flux 
emitted by the slab can be calculated. See CEP06 for 
more details and the derivations of the original plane- 
parallel equations to which we make reference below. 

2.1. Theoretical/ Analytical Expressions for the Net 
Radiative Bracket 

In CEP06 equation 7, Elitzur & Asensio Ramos de- 
rive a purely analytical expression for the net radiative 
bracket in a plane parallel slab, based on the formal solu- 
tion of the radiative transfer equation and the definition 
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of the net radiative bracket: 



p(r) 



1 - 



S(t)dt / <5> 2 {x)dx / e -|T-*l*(«)//i 
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where r and t refer to optical depths for a specific 
wavenumber, with t% being the overall optical thickness 
of the slab, S(t) or S(t) to the source function for a given 

line, {S Vo , — -5 — A2in J} ) $(x) is a dimensionless line 

profile, and fi — cos0, where 9 is the angle of a given ray 
measured from the normal to the plane. 

We use a spherical analog to this theoretical expres- 
sion (i.e. as opposed to the discretized expression they 
introduce later involving a number of "zones") for the 
net radiative bracket: 

p(r(r,M)) = 



1 - 



4tt 
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T(r,e,4>,Q) 



2S(T(r,6,4>,n)) J 
(x) ^ e -|r(r,fl,*,n)-t(r,M,n)l*(»)) dx dt 



S(t(r,M,fi) 

dn (3) 



where p(r(r, 9, <j))) refers to the net radiative bracket at 
any point labeled by the coordinates (r, 9,<f>) in spherical 
coordinates with the origin at the center of some sphere 
of interest. (In our study, the sphere is centered on the 
comet nucleus; it could be the center of any spherical 
astronomical object for any given case.) 

As in CEP06, = (tt A^)" 1 / 2 e'* 2 is a dimen- 

sionless line profile, normalized so that J $>{x)dx = 1, 
where x — — i;q) / Avp is the dimensionless line shift 
from line center vq and Avjj is the doppler line width, 
Vq/ c{2kT /m) 1 / 2 . In the present work, we neglect doppler 
offsets between different locations in the coma with rel- 
ative velocities. This is due to our work being primarily 
motivated by the Deep Impact observations of the near- 
nucleus coma (see Feaga, et al. 2007), for which gas 
tends to be released, broadly speaking, towards one di- 
rection/side of the nucleus. For this situation, doppler 
shifts are likely to be less than line broadening. In fu- 
ture work, we plan to include a more generally accurate 
treatment. 

Both T(r,6,(j),dfl) and t(r,9,(f>,dSV) refer to optical 
depths and S(r(r, 6, <f>, dCl)) to the source function as 
viewed from the coordinates (r, 8, cf>) along a "pencil" of 
solid angle dQ. (See Fig. |J) 

The optical depth along a "pencil" (or cone) of solid 
angle d£l is, of course, highly dependent on the particu- 
lar direction of a given solid angle element d£l, i.e. de- 
pendent on 9',4>', the direction of a vector centered on 
(r, 9, 4>) pointing along the centerline of dfl. (See Fig. [2j) 

2.2. Discrete Expressions for p 

The above CEP06 expression for p(r) (Equation [2| is 
turned into a discrete expression by dividing the slab into 
multiple "zones" (labeled with indices i This yields 




Fig. 2. — Sphere showing a point (r, 0, 
elements dQ± and dfl2- (Viewed from the 

CEP06 equation 14: 



)) and two solid angle 
-y direction.) 
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This expression for p can be calculated for each zone 
(and each wavenumber/frequency) and the resultant p's 
included in the equations of statistical equilibrium, which 
are then solved. (See also CEP06 equations 32-36, and 
accompanying CEP06 text, for a more complete discus- 
sion.) 

Similar to CEP06 Equation 14, for the purposes of our 
integration of a discretized p l the sphere is divided into 
spherical shells (analogous to the plane-parallel zones of 
the original CEP) where i (or j) is a shell index. (See Fig. 



The integration is broken down into a sum of integrals 
ong different "cones" of solid angle (the aforementioned 
dfi's). Note that although d£l may seem somewhat con- 
ceptually analogous to d\i in the plane-parallel case, we 
cannot use the convenient exponential integrals as in the 
plane-parallel situation, but instead must integrate along 
each d£l separately. 



5> J / dr 



dt 



\x)dx 



-|r(n)-t(n)|*(x) 



(n) 

dn 



(5) 

Here we have dropped the coordinate subscripts r,9,(f> 
for clarity /simplicity and use shell indices instead (where 
some shell i contains the point defined by r,9,tfi). Note 
that the "z" in the summation, the maximum number of 
shells along a given d£l, will be different for each dQ. For 
each "cone" of solid angle we will sum S 1 - 5 (f2)e - ' T-t ' from 
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Fig. 3. — Example sphere showing division into 4 shells, i = 1..4 
and solid angle dQi for which 2 = 2 with indices j = 1,2, as per 
Equation [5] (Viewed from the — y direction.) 



"z" (the outermost shell) to i - 
adjacent to shell i. 



1, where i + 1 is the shell 



2.3. Further steps towards implementation 

In actual practice (i.e. computer implementation), the 
integration of the source function of dQ over 4n steradi- 
ans will also be done by a discrete summation. 

Along any particular element of dfl originating in shell 
i, this sum will be: 



^^'(d^Xl-e- 

Z 



(do.)^ e -T 3 - 1 - i (dn) 



) n 



(dfi) 



(6) 



Note that, unlike the plane-parallel case, each dr must 
be calculated explicitly from the geometry and local 
molecular energy level populations for each shell along 
the "pencil": 



dr 3 {dVL) = a j (dn) x ds j (dn), 



(7) 



where a 3 (df2) is the absorption coefficient in region j and 
ds J (dfl) is the distance through that region. Note that 
this distance may vary over the width of a given dQ but 
we can either approximate using the centerline of dfl or 
derive a ^'-dependent expression and calculate a proper 
integration to get a mean value over the region. (In our 
implementation we use the former option, for the sake of 
simplicity.) 

This integration over solid angle is more tedious than 
in the plane-parallel case (which was able to use exponen- 
tial integrals over a zone's dr) and more computationally 
costly, but straightforward enough to be feasible. 

Turning this all into a fully discrete expression for p l 
(for shell "z") we get 



P = 
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Where Z is the maximum numbered spherical shell and 
all the quantities indicated by the subscript uj are depen- 
dent on a particular direction viewed from a given point 
(or shell i, in a spherically symmetric case). 

2.4. Asymmetric case: Incident Radiation 

Our motivating interest in this study is comets' comae, 
which are not spherically symmetric cases. 

To adapt CEP to this asymmetry we divide each shell 
further into "regions." In the case of comets, one source 
of asymmetry is incident solar radiation coming from one 
direction (outside the outermost shell). Therefore, the 
natural way to further divide shells into regions is along 
lines parallel to the direction of solar radiation (i.e. the 
center-of-comet-to-center-of-Sun line, which we will ar- 
bitrarily label as the z-axis.) Thus we superimpose a set 
of co-axial cylinders (with radii equal to corresponding 
shell radii, for the sake of simplicity) on the shells to 
divide the coma into regions bounded by two cylinders 
and two spherical shells. (Note that some regions, specif- 
ically those along the z — plane perpendicular to the 
solar radiation, are only bounded by an inner cylinder 
and outer sphere. See Figs. |i]&[5j) 




Fig. 4. — 2-D cross section of sphere on y = plane, viewed from 
the — y direction, showing division into 4 shells, with superimposed 
cylinders along solar (2) direction. 

These regions form annuli or rings of unusual, but eas- 
ily envisioned, cross-sections. (See Fig. [5]) 

Incident solar radiation is parallel to the z-axis (due to 
our choice of the z direction). Hence each ray of sunlight 
travels along the axial direction of a specific cylinder. For 
each region, the solar radiation absorbed is calculated 
based on the relevant optical depths along that direction 
(the dr's) of those regions in the same cylinder that are 
closer to the Sun than the given region. This is similar 
to the case of external radiation described in CEP06 Ap- 
pendix A, equations A3 and A4, in which we simply set 
fio = 1) due to the above constraint of cylinders being 
co-axial with the incident solar radiation: 
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Fig. 5. — Two examples of different possible shapes of 3D annuli 
formed by intersecting spheres and cylinders. 

J| = Je^bir*)- 7(^-1)] (9) 
where 

where J* is the average over a region i of J\, the con- 
tribution of external radiation, J e , to the mean intensity 
of the region, which is to be included in the equations 
of statistical equilibrium (Equation [TJ by addition to the 
"B" term. 

From a purely radiative standpoint, assuming that 
within each region there exists uniform density, temper- 
ature and other physical conditions, the radiative ex- 
citation of molecules (hence, the emission and absorp- 
tion) in each region/annulus should be equal through- 
out the region. In our expanded CEP implementation, 
these regions are analogous to zones in the plane-parallel 
CEP. Each region's radiative effect or contribution to 
each other region, i.e. the net radiative bracket, must 
be calculated. Note that self-irradiation from around an 
annulus must also be taken into account, as well as irra- 
diation from other regions. Once this calculation is done, 
the entire region will be equal with respect to radiative 
processes (i.e. there is symmetry around the z-axis). 

2.5. Further Asymmetry: Coma Morphology 

If models of distantly observed comets were all we 
needed, this might be sufficient. But we are motivated by 
the desire to better understand Deep Impact & EPOXI 
spectral observations that have very high spatial resolu- 
tion around the comets' nuclei. (See e.g. Feaga et al. 
2007 and Feaga, et al. 2011.) Therefore the above ra- 
diative treatment alone is insufficiently asymmetric to 
fully model a cometary coma, when coma morphology 
is included in the model. The inclusion of morphology 
undoes the aforementioned symmetry around the z-axis 
within each annulus/region. These observations are one 
of the primary motives for this study, and therefore these 
morphological asymmetries must also be dealt with ap- 
propriately in this model. 

To model morphological features, we use a cone shape 
superimposed over the aforementioned divisions into re- 
gions. (Other geometric shapes could also have been 
used. We chose to implement a cone due to its sim- 
ilarity in shape to many observed coma features.) A 
cone of arbitrary orientation and size with its vertex at 
the center of the sphere creates intersections with the 
above-described regions. Each of these is then added as 
a sub-region, which may have different properties from 



the surrounding (or subsumed/replaced) region. 

Each sub-region can posess different initial conditions 
from the surrounding region. Thus morphological fea- 
tures, which by their nature tend not to be axisymmet- 
ric around our z-axis, can be included in the model. It 
should be noted that these sub-regions are only included 
as necessary. Thus for those annuli that do have constant 
axisymmetric conditions (i.e. no interesting morpholog- 
ical features impinging on them), we can save time and 
memory computationally by leaving them undivided, as 
they would have been originally. 

2.6. Implementation: Our Algorithm Described 

Given the above geometric divisions, we have imple- 
mented Asymmetric Spherical CEP as follows. For each 
region (or sub-region, as applicable) we take representa- 
tive population values from the coma integration and use 
these as the basis of an "initial guess." We then make 
an immediate improvement to the initial guess values 
by recalculating each region's populations (individually) 
taking into account the attenuation of incident solar radi- 
ation by intervening regions in the solar direction (as per 
Equation [9j|. These recalculated populations are the val- 
ues we then use as the initial guess (required by the im- 
plementation of Newton's method in Press, et al. 1992) 
for CEP calculations. Based on these populations we cal- 
culate the necessary source functions, delta-taus, and net 
radiative brackets, "p," as above, for each wavenumber 
(or line, transition, etc.) Following the above discretized 
equation Q for each region, which in this context wc 
will call the "recipient" region, we iterate over all other 
regions to calculate their contribution to the recipient's 
p. Each other region's contribution is essentially its own 
source function attenuated over the optical depth of all 
intervening regions along the line of sight between itself 
and the recipient region, integrated over (or, to simplify, 
multiplied by) the solid angle subtended by one region 
from the other. This is then divided by the recipient re- 
gion's source function and optical depth (along the given 
line of sight). 

To implement this in a practical algorithm of manage- 
able complexity, we make several simplifying approxima- 
tions. 

Due to the z-axis symmetry that exists (before 
adding sub- regions), we can partially simplify to a two- 
dimensional diagram in which a region is represented by 
the cross-section of the annulus in the (arbitrarily cho- 
sen) y = plane. We calculate a region's "centroid" , i.e. 
the centroid of its 2-D projected area in this plane which 
(in our approximation) corresponds to a point (r, #, <f>) in 
spherical coordinates. 

For (cone shaped) sub-regions, which in general do not 
have their centerline on the X-Z plane, we must use a 
different centroid. We use the midpoint along the cone's 
centerline (within region boundaries). 

We also choose a series of points distributed evenly 
along circles parallel to the X-Y plane around each re- 
gion, which will be the "starting points" for calculation 
of lines of sight (which will terminate at the centroids). 
These are chosen by rotating a region's centroid around 
the z-axis by multiples of some angle that depends on the 
size (radius) of the region. The choice of angle is such 
that the larger a region's size, the more starting points it 
will have, and thus the region will be divided into more 
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elements of solid angle. 

We use the line of sight between the "centroids" of re- 
gions and this series of starting points to calculate the 
contributions of every other region (or the region to it- 
self) to a given recipient region's p. We calculate the op- 
tical depth of each intervening region, along the line of 
sight, based on the molecular population levels of the in- 
tervening regions. (See Fig. [6]) These "integration lines" 
encapsulate the main part (within the square brackets in 
Equation Is!) of the calculation of p. 





Fig. 6. — A two-dimensional view from above (i.e. the +z di- 
rection) of examples of integration lines in the X — Y plane. Four 
lines originating in the i = 1 region and ending at the centroid 
of the i = 2 region are shown (in red in the online version), with 
corresponding division of the i = 1 circle (shown by dashed grey 
lines). Eight lines originating in the i = 4 region and ending at the 
centroid of the i = 3 region are shown (in blue in the online ver- 
sion), with corresponding division of the i = 4 annulus (also shown 
by dashed grey lines). One example of a corresponding df! is also 
shown with dotted-dashed lines (in green in the online version). 
Note that this 2-D diagram only shows horizontal cross-sections of 
regions, and so regions and shells/annuli are essentially indistin- 
gushable in this diagram. 

We approximate the solid angle subtended by the inte- 
gration lines from another region's centroid by integrat- 
ing dQ = d<j> svaBdO from zero up to the mean value of 
the angles between the starting point of that region, the 
centroid of the recipient region and the multiple "cor- 
ners" of that region around the starting point. (See Fig. 
[7 ) Effectively, this gives a solid angle between regions 
of 2tt(1- cos0 mea „). 

Note that due to these approximations, the sum of solid 
angles over all integration lines between a region and all 
regions in a given shell exterior to the region is not neces- 
sarily constrained to exactly equal 4n, as it should be in 
reality. Therefore, in calculating a "p" value, we sum the 
solid angles involved and average over that sum instead 
of 47r steradians (as Eqns. [3] & p] dictate we should do). 

In the limit of arbitrarily small (and numerous) re- 
gions, these approximations would approach a physical 
situation of arbitrarily precise accuracy. Thus we main- 
tain the "exactness" of the CEP method. 

Unlike the plane-parallel situation, the flux exiting the 
surface of the coma (or other sphere of interest) is not 




Fig. 7. — Example illustrating calculation of mean angle. Lines 
originating from an angular slice of a region's "corners" and ter- 
minating at another region's centroid are shown in black. The line 
between the centroid of the "recipient" region and a "start point" 
of the other region (the largest dot), corresponds to the relevant 
integration line (and is shown in red in the online version). The 
integration line and each of the other eight lines define the angles 
that are averaged together to get the mean angle <9 me an used to cal- 
culate the solid angle subtended by one region as viewed from the 
other's centroid. Note that not all regions will have eight "corner 
points." 

simply a single value (per wavenumber) that has been 
integrated over angle. In the spherical situation, the re- 
sultant intensities form a two-dimensional mapping (in a 
planeperpendicular to the observer's line of sight. See 
Fig. Hf) 




Fig. 8. — As per Fig. [4] above, with observer plane also shown, 

rk-. 



aligned perpendicular to 
online version.) 



-axis. (A color image is included in the 



In our implementation, this plane is specified by rota- 
tion angles, 0, (f> and ip with the comet's center at the 
origin, and is assumed to be at a distance > R coma , the 
maximum radius of the comet's coma. We can also spec- 
ify the density of and interval between points on this 
plane for which the output intensities will be calculated. 
Each point in this planar mapping shows the intensity 
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(or surface brightness) integrated along a specific line of 
sight, perpendicular to the plane, through the coma from 
one side to the other (again, for each wavenumber). 



Isurf 



J2 Si^ii 1 



i=l 



3=1 

' n • 

j=i— 1 



(10) 



Where Si is the source function of a region i, Aui is the 
line width of wavenumber/frequency v in region i, and 
t % , or r J , represents the optical depth of wavenumber v 
in region i or j along the relevant line of sight. Indices 
i and j run from 1 to 2, where z equals the number of 
regions along a given line of sight. 

Thus the spherical CEP algorithm produces results 
that could be described as a four-dimensional data "hy- 
percube" : for each point in the above 2-D mapping, there 
is a complete (2-D flux vs. wavenumber) spectrum. This 
data can then be presented in multiple formats. Several 
forms of data presentation for simulating observations 
are described in the following section. 

This is also precisely the output needed to compare 
with the Deep Impact and EPOXI observations that have 
been displayed as two-dimensional brightness maps for 
specific wavelengths or bands. 

3. SOME PRELIMINARY RESULTS: OBSERVABLES FOR 
DISTANT COMETS 

We present here examples of model results for three dif- 
ferent production rates of carbon monoxide which could 
be potentially useful for distant (e.g. ground based or 
orbital telescope) observers of comets. These are mod- 
eled using a spherical coma, without any morphological 
features but including optical depth effects both with re- 
spect to incident solar radiation within the coma and 
with respect to emergent "observed" radiation. 

The output data from the CEP model can be presented 
in various ways. Here, we show an example of a band to- 
tal brightness map (analogous to Feaga, et al. 2007, but 
for the entire coma), radial profiles of brightness, col- 
umn density, and g-factors for various azimuthal angles. 
(These could also be done for individual spectral lines, 
but in the interests of space and avoiding complexity we 
have not presented such results here.) We also present 
spectra integrated over different "aperture" sizes. The 
band total brightness is more likely to be similar to ac- 
tual observations, but high resolution spectra are possi- 
ble, even from ground based telescopes (see e.g. DiSanti, 
et al. 1999 and DiSanti, et al. 2001), in particular for 
comets close to Earth, which might more closely resemble 
the latter form of model results. 

We also demonstrate the potential usefulness to ob- 
servers of the ratio of the total brightness of the P branch 
to that of the R branch of CO to determine whether ob- 
servations include significant optical depth effects. This 
may be measurable even with relatively poor spectral 
resolution. 

Many of the model "input values" (see Table [lj have 
been chosen so as to facilitate comparisons (of our op- 
tically thin cases) with other earlier models. Model pa- 
rameters that are the same for all the following exam- 
ples are: Solar distance = 1 AU. Solar flux (over the 
CO band) is 2.5 x 10 13 photons cm~ 2 s _1 (cm -1 )" 1 , as 



per Labs & Neckel 1968 (and as used by Chin & Weaver 
1984 and Weaver & Mumma 1984). Gas expansion speed 
is a constant 0.8 km s _1 and the initial gas temperature 
at the surface is 200 K. Qh 2 o = 10 x Qco, as in Chin 
& Weaver 1984. As mentioned above, the radial tem- 
perature profile closely follows Combi's 1989 model (see 
Fig. [IJ but scaled to the initial gas temperature at the 
surface, T slir y ace . 

The coefficients for CO-H 2 collisions (assumed to be 
the dominant source of collisional excitation) are as per 
Chin & Weaver 1984: only rotational excitation and de- 
excitation are considered. (Vibrational cross sections 
are about 5 orders of magnitude smaller. See Weaver 
& Mumma, 1984, Table 2.) C = n H20 crv, where v is 
the average relative speed of the molecules (cm s _1 ), 
n H20 is the number density of H2O (cm -3 ) and a is 
the collisional cross section of a given transition of CO 
(cm 2 ). The last value is based on a total cross section of 
a tot — 1-32 x 10~ 14 cm 2 , which is apportioned between 
AJ's up to 6 as per Chin & Weaver's Table 1, which we 
reproduce here in our Table [2j 

3.1. Brightness Maps 

We present here one example of a brightness map of a 
modeled coma (see Fig. [9]) . This format of output pre- 
sentation is most similar to the radiance maps of Feaga, 
et al. 2007 and A'Hearn, et al. 2011. For a spherical 
coma with no morphological features, it is rather unin- 
teresting. It is nevertheless included here as a demon- 
station and used to illustrate the azimuthal angles of 
the radial profiles presented below with the addition of 



28 



case, 



overlaid lines. Note that for the Qco = 10 
there is some difference in brightness noticeable to the 
eye between the sunward side (azimuthal angles nearer 
to zero) and the anti-sunward side, especially in the near- 
nucleus portion of the image. This is due to the optical 
depth along the solar direction reducing the excitation 
and emission from one side of the coma to the other. 

3.2. Radial Profiles: Brightness, Column Density and 

g-factors 

Abundances of cometary species are frequently calcu- 
lated from observed fluxes using fluorescence efficiencies, 
or g-factors. In an optically thin case, the brightness of a 
given line or band is proportional to the column density 
of the relevant molecule. In such cases Fb an d = 9band X N 
and g band = J2band A u x n « where F band is the band to- 
tal flux (or brightness), gband the band g- factor, N the 
total column density, A u the Einstein A coefficient for 
the relevant transition originating in upper level "u" and 
n u is the column density of the population of a specific 
upper level "u" (which in our model is numerically ap- 
proximated as the sum over all regions along a line of 
sight of the fractional population of level "u" times each 
region's column density). 

However, large optical depths will spoil this simple 
linear relation between column density and brightness. 
With radiative transfer modeling, it is possible to get a 
calculated g-factor {g band = J2band A u x n u ) from the 
model and the "effective g-factor", g e ff = Fb an d/N, 
which is the actual ratio of brightness to column den- 
sity. 

In Figs. 10 |TTJ and 12 we present all these val- 
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TABLE 1 

Model input parameters for models of theoretical example comets for CO. 



Mean Nucleus Radius 


3 km 


sur f ace 


200 K 


Expansion speed V e xp 


0.8 km s" 1 


Qco 


io^.io^.io 38 *- 1 


Band Center Wavenumber 


2149 cm,- 1 


Band Center Einstein A 


33 s" 1 


Highest Rot. Level, Jmax 


20 


Solar flux 


2.5 X 10 13 photons cm"' 1 s -1 (cm -1 ) -1 


0~rot 


1.32 X 10" 14 cm 2 



TABLE 2 

Reproduction of Chin & Weaver's Table of CO-H2O 
Collision Cross Section Information. 



A J — Supper 


- Slower Fraction of Total De-Excitation 


1 


0.34 


2 


0.25 


3 


0.20 


4 


0.10 


5 


0.07 


6 


0.05 


>6 






Total cross section is always o> ot = 1.32 X 10 ci 
Excitation is derived from de-excitation using 
detailed balance. 

CO 1-0 Bond Tot. Brightness Mop, 0=1 0™s~' (Phose angle = 90°) 
600 ■ 




200 400 600 



ergs s" cm" sr" 



Fig. 9. — Example of a band total brightness map for the CO 1-0 
band, for the inner ±600 km near the nucleus of a theoretical comet 
with Qco = 10 28 s — , viewed from phase angle = 90° . Overlaid 
radial lines indicate the orientation of azimuthal angles in radial 
profiles below. The sunward direction is up, i.e. azimuthal angle 
= 0° . Note that within the ±600 km field of view, the brightness 
never reaches zero (even where it appears to be totally dark). This 
image is in color in the online version of the article and the colors 
of the azimuthal angles correspond to those in subsequent radial 
profiles. 

ues together as radial profiles, for theoretical comets 
of three different production rates, Qco — 10 26 , 10 27 , 
and 10 28 s _1 . (All observed at 1 AU, at a phase an- 
gle of 90° and multiple azimuthal angles. In all cases 
Qh 2 o = 10 X Qco-) 

In our results, we typically see that <? e // does tend 
towards the calculated g-factor values at larger impact 
parameters, where optical depth effects are negligible, 
as would be expected. (The actual numerical values of 
the "asymptotic" band g-factors produced by our model, 
2.4 x 10 _4 s _1 per molecule for CO at 1 AU, also agree 
well with other published values such as those calculated 



by Chin & Weaver 1984, Crovisier & Le Bourlot 1983, 
and Weaver & Mumma 1984.) The actual radii at which 
this convergence occurs depends primarily on the produc- 
tion rate of a comet. We can use the distance at which 
9eff = Q-9gthin as a very rough measure of the point 
where a coma can be considered to transition from op- 
tically thick to thin. For the "thin" and "intermediate" 



coma models (Qco = 10 2B and Qco = 10 27 s_1 ) the con- 
vergence is fairly close to the nucleus, within ^100-200 
km. But for the "thick" model (Qco = 10 28 s _1 ) with its 
high production rate, the "optically thick regime" can ex- 
tend as far as O(10 3 ) km, which can be spatially resolved 
even in some remote observations. Note that at radial 
distances very near the nucleus, worrying about optical 
depth effects may be relevant even for lower production 
rates. 



3.3. Radial Profiles: 



Phase and Azimuthal Angular 
Variations 



Another optical depth effect is variation of brightness 
(and corresponding g-factor) with varying angles, both 
phase angle (of observer) and azimuthal angle within any 
given observation. 

The radial profiles in Figs. 



10 12] demonstrate the az- 



imuthal variation for a single phase (observing) angle. 
Not surprisingly, the farther a radial profile line is from 
from the sunward direction, the less bright it is and the 
lower its g-factors for a given distance from the nucleus. 
However, the degree to which the actual g-factor varies 
is of note. Even for the only moderately thick case of 
Q co = 10 27 s~ 1 for radii < 100 km, there is a difference 
of as much as ~10% between sunward and anti-sunward 
directions. The effect is even more pronounced for the 



thicker cas e of Qc o 
In Figs 



10 



2H C 



12 



12||16 we present profiles of brightness (and 
column density) tor model results observed from different 
phase angles, ranging from 0° to 180° , for the optically 
thicker case of Qco — 10 28 s _1 . Each plot includes mul- 
tiple azimuthal angles, which would all be visible simul- 
taneously in a wide field observation (i.e. including the 
entire coma) from each given phase angle. (The ab ove 
plot for 90° phase angle for Q C o = 10 28 s~\ Fig 
should be considered part of this series as well.) 

Observations with a slit spectrometer with sufficient 
spatial resolution (see, e.g. DiSanti, et al. 1999 ) might 
observe along one specific azimuthal angle and thus see 
possible variations along the slit with sufficiently high 
spatial resolution. 

The most obvious effect seen at a glance in these figures 
is the spread among azimuthal angles for a given phase 
angle. As would be expected, the 0° and 180° phase an- 
gles (sunward and anti-sunward) have no real azimuthal 
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CO 1-0 Band Tot. Q c0 = 1x1 0"s"': Brightness, Col. Density 
& Ratio vs. R (Phase angle = 90°, Multiple Azimuthal angles) 



CO 1-0 Band Tot. Q c0 = 1x1 ! V: Brightness, Col. Density 
& Ratio vs. R (Phase angle=90°, Multiple Azimuthal angles) 



1.0x10" 
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4x10" 
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Az. = -135 
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10 3 
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Fig. 10. — For Qco = 10 26 s _1 . Upper frame: Radial profile of 
band total Brightness vs. R (impact paramater) for Phase angle 
= 90° and for multiple Azimuthal angles. Profiles of azimuthal 
angles (indicated by color coding in the online version) show no 
variation for this case and overlap, appearing indistinguishable. 
Column density is included as the bold solid line (red in the online 
version) using a different y-axis scale. Lower frame: g-factors. 
Both the calculated g- factor (the higher group of lines) and g e ff, 
the observed brightness over column density, plotted with matching 
styles (colors in the online version) for each azimuthal angle (which 
also match those in the upper frame). Profiles for 0° , ±45° , and 
90° overlap each other almost entirely. 



variation. From phase angle 45° to 90° to 135° there is a 
progression: the azimuthal lines get spread out farther 
from each other, as well as noticably dropping in bright- 
ness for those lines farther from the sunward side. 

With respect to total brightness, the phase angles 
0° , 45° and 90° , produce roughly equal peak bright- 
ness for their strongest azimuthal profiles (the more sun- 
ward directions, at the nucleus grazing radius) of about 
~ 8 x 10 11 photons s -1 cm 2 sr _1 For the 135° phase 
angle, the peak values are about ~ 6 x 10 11 . Most sig- 
nificantly, for the 180° phase angle, the peak values are 
about ~ 3 x 10 11 - only about half as bright as at other 
phase angles. 

With respect to "effective" (or "observed") g-factors, 
the minimum values in the most optically thick regions 
(again, at the nucleus-grazing radii) also show a trend 
from sunward to anti-sunward. For the 0° , 45° , and 
90° phase angles, the minimum values of the ratio of 
the band total flux over the column density are roughly 
1.3 x 10 -5 photons s _1 molecule -1 . From 90° , through 
135° and down to 180° the values go down monotonically 



10 3 
Radius (km) 

Fig. 11. — For Qco = 10 27 s _1 . Upper frame: Radial profile of 
band total Brightness vs. R (impact paramater) for Phase angle 
= 90° and for multiple Azimuthal angles. Profiles of azimuthal an- 
gles (indicated by color coding in the online version) for 0° , ±45° , 
and 90° overlap each other almost entirely and other brightness 
profiles are almost indistinguishable. Column density is included 
as the bold solid line (red in the online version) using a different 
y-axis scale. Lower frame: g-factors. Both the calculated g-factor 
(the higher group of lines) and g e //, the observed brightness over 
column density, plotted with matching styles (colors in the online 
version) for each azimuthal angle (which also match those in the up- 
per frame). Profiles of azimuthal angles 0° , ±45° , and 90° overlap 
each other almost entirely and are almost indistinguishable. 



to a minimum of about 4.5 x 1CP 6 . This trend is due 
to a combination of two optical depth effects. The first 
is attenuation of incident solar light from the sunward 
to anti-sunward sides of the coma, leading to less flu- 
orescent pumping, and thus less emission, towards the 
anti-sunward direction. Second, whatever emission there 
is is more likely to "escape" the coma over shorter opti- 
cal depths - i.e. closer to where it is emitted. Thus the 
already greater emission of the sunward regions is also 
more likely to be observed along azimuthal directions 
closer to sunward. 

However, the similarly located values for "actual" cal- 
culated g-factors do not follow a similar simple monotonic 
trend. The greatest values for a given phase angle rise 
from 0° through 45° and peak for phase angle 90° . From 
90° , through 135° down to 180° they fall through the 
same values, crea ting a symmetric peak around 90° . For 
90° (see Fig. " 



12 1 there are two "clusters" of lines. The 
higher one corresponding to more sunward azimuthal an- 
gles (between ±45) and the lower one to other angles. 
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CO 1-0 Band Tot. Q c0 = 1x1 0'V: Brightness, Col. Density 
& Ratio vs. R (Phase angle = 90°, Multiple Azimuthal angles) 




10' 10 2 10 3 10* 10 5 

Radius (km) 



Fig. 12. — For Qco = 10 28 s _1 . Upper frame: Radial profile of 
band total Brightness vs. R (impact paramater) for Phase angle = 
90° and for multiple Azimuthal angles. Profiles of azimuthal angles 
(indicated by color coding in the online version) for 0° , ±45° , and 
90° overlap each other almost entirely. Column density is included 
as the bold solid line (red in the online version) using a different 
y-axis scale. Lower frame: g-factors. Both the calculated g-factor 
(the higher group of lines) and g e //, the observed brightness over 
column density, plotted with matching styles (colors in the online 
version) for each azimuthal angle (which also match those in the 
upper frame). 

The upper cluster's (minimum) values are actually the 
highest values in this comparison, about ~ 6.1 x 10~ 5 . 
The lower cluster of values at 90° is still greater than 
either angle 0° or 180° , which are the lowest for any 
phase angle, notwithstanding the azimuthal spread for 
the other phase angles. (The calculated values also show 
considerably more spread among azimuthal angles than 
the profiles of flux over column density.) This symmet- 
rical and non-monotonic pattern is less intuitive than 
the trend of <? e // above. Yet it is clearly understood in 
light of the fact that these values are based only on the 
actual population distributions in different regions and 
do not include optical depth effects on the emergent ra- 
diation. Thus observing from azimuthal angles 0° and 
180° are sampling exactly the same lines of sight and re- 
gions' populations, including the darkest (i.e. least ex- 
cited populations) of the anti-sunward side of the coma. 
The same is essentially true for 45° and 135° (due to sym- 
metry around the z-axis) but they do not sample the 
darkest parts of the anti-sunward directions (and the 
differences in populations are less between outermost re- 
gions on the sunward side among azimuthal angles be- 



CO 1-0 Band Tot. Q c0 = 1x1 0'V: Brightness, Col. Density 
& Ratio vs. R (Phase angle = 0°, Multiple Azimuthal angles) 



- 






Az.=0 (Sunward) 




A?. =45 




.... Az. = 90 




Az.= 135 




A7 =-1Rn fAnti-^unwnrr-n 




A2.--135 




-Az. = -90 




.... Az. = -45 







6x10 ie 



10' 10 2 10 3 10" 10 5 

Radius (km) 



g-factors 





2.41x10"* 
2.17xl0" 4 










lotons s~' molecule"') 




Colculoted (=E A„x 


N„) 


A>. = 45 
.... Az. = 90 
Az.= 135 

A7 =-1Rn (Anti-Sun 

Az. = -135 

j<\z. = -90 

Az. = -45 


Q. 












o 












u 
o 












cn 


4.46x10" 5 


^- 










1.27X10" 5 




--- --""Observed" g, ( i( = B 







10' 10 2 10 3 10" 10 5 

Radius (km) 



Fig. 13. — For Phase angle = 0° , Qco = 10 28 s _1 . Upper frame: 
Radial profile of band total Brightness vs. R (impact paramater) 
for Phase angle = 90° for multiple Azimuthal angles. Profiles of 
azimuthal angles (indicated by color coding in the online version) 
show no variation for this case (as should be expected) and over- 
lap, appearing indistinguishable. Column density is included as 
the bold solid line (red in the online version) using a different y- 
axis scale. Lower frame: g-factors. Both the calculated g-factor 
(the higher group of lines) and 9 e //, the observed brightness over 
column density, plotted with matching styles (colors in the online 
version) for each azimuthal angle (which also match those in the 
upper frame). 

tween ±45 - they are all experiencing direct solar illu- 
mination). For 90° the higher cluster is sampling from 
more excited and higher emmision populations than the 
lower cluster, and consistently so all the way along their 
lines of sight (which is not true for 45° and 135° ). Thus 
the sunward cluster of lines for 90° is the brightest seen, 
and the anti-sunward cluster values falls between values 
of the 0° or 180° and the 45° or 135° lines. 

Lastly, the profiles for Qco — 10 28 s _1 have a "bump" 
in brightness in the vicinity of ~ 100 ~ 1000 km, most 
easily visible for the 180° plot, but also present for the 
other phase angles (and growing in size from 0° up to 
180° ). This is mostly due to the temperature profile 
reaching its minimal values at these radii in conjunction 
with the higher density of the Qco — 10 28 s _1 case. The 
higher density leads to this still being a collsionally dom- 
inated regime, and the low temperatures lead to the low- 
est population levels being most highly populated. These 
levels also have the highest Einstein A values, thus lead- 
ing to higher overall number of photons emitted for the 
same number of molecules. The effect is greatest for the 
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CO 1-0 Band Tot. Q c0 = 1x1 2 V: Brightness, Col. Density 
& Ratio vs. R (Phase angle = 45°, Multiple Azimuthal angles) 
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Fig. 14. — For Phase angle = 45°, Qco = 10 28 s _1 . Upper 
frame: Radial profile of band total Brightness vs. R (impact para- 
mater) for Phase angle = 90° for multiple Azimuthal angles. Pro- 
files of azimuthal angles (indicated by color coding in the online 
version) arc slightly displaced for viewing purposes. Column den- 
sity is included as the bold solid line (red in the online version) 
using a different y-axis scale. Lower frame: g-factors. Both the 
calculated g-factor (the higher group of lines) and g e //, the ob- 
served brightness over column density, plotted with matching styles 
(colors in the online version) for each azimuthal angle (which also 
match those in the upper frame). 

180° view due to a cumulative effect - the lines of sight 
all sample the most dense and cold regions at these radii. 
The same extreme effect is not seen for the 0° phase an- 
gle due to the overall greater fluorescent excitation of the 
sunward regions dominating it. (Note that the difference 
in total brightness between the two azimuthal angles at 
these radii is about a factor of two.) 

3.4. Aperture Averaged Spectra 

If one is observing with high spectral resolution but 
low spatial resolution, the spectra observed will be the 
sum of as much of the coma as fills the field of view. 
To model this, we have simulated "aperture averaged" 
spectra, where the "aperture" controls the area of the 
coma sampled. Our apertures are square boxes and are 
all centered exactly on the center of the comet, and sam- 
ple a nucleus centered area equal to the square of the 
"aperture size" over which we average the brightness. 
We present a series of apertures from 2 x 10 1 km (very 
near the nucleus) through 2 x 10 5 km (the whole coma) 
for each of the three production rates. (All these exam- 
ple spectra are modeled at a heliocentric distance of 1 



CO 1-0 Band Tot. Q c0 = 1x1 ! V: Brightness, Col. Density 
& Ratio vs. R (Phase angle=135°, Multiple Azimuthal angles) 
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Fig. 15.— For Phase angle = 135°, Qco = lO 28 * -1 . Upper 
frame: Radial profile of band total Brightness vs. R (impact para- 
mater) for Phase angle = 135° for multiple Azimuthal angles. Pro- 
files of azimuthal angles (indicated by color coding in the online 
version) arc slightly displaced for viewing purposes. Column den- 
sity is included as the bold solid line (red in the online version) 
using a different y-axis scale. Lower frame: g-factors. Both the 
calculated g-factor (the higher group of lines) and g e ff, the ob- 
served brightness over column density, plotted with matching styles 
(colors in the online version) for each azimuthal angle (which also 
match those in the upper frame). 

AU and phase angle 90° . Qh 2 o = 10 x Qco, as in Chin 
& Weaver 1984.) 

Band shape for apertures including the outer coma (ap- 
proximately 10 4 < R ap < 10 5 km) does not change signif- 
icantly for different production rates. The total bright- 
ness for this regime increases approximately in linear pro- 
portion to production rate. This is due to spectra with 
such large aperture sizes being dominated by the fluores- 
cence dominated optically thin outer coma with optical 
depth effects playing a minimal role. (But not entirely 
non-existent: note the small, < 6% reduction in g-factor 
with higher production rate for R ap = 2 x 10 4 km.) 

In the inner coma, however, optical depth effects can 
be very striking. The "thickest" spectra (for higher pro- 
duction rates and smaller R ap ) have remarkably altered 
band shapes from the optically thin spectra. 

First, however, a word about changes that are not 
specifically caused by optical depth effects. It is clear 
that there is much variation, even within a given pro- 
duction rate, from large to small aperture sizes. Not all 
of this is due to optical depth effects. Even in our op- 
tically thin case, the band shape changes noticably in 
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CO 1-0 Band Tot. Q c0 = 1x10 ! V 
& Ratio vs. R (Phase angle=180°, 



Brightness; Col. Density 
Multiple Azimuthal angles) 
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Fig. 16.— For Phase angle = 180°, Q C o = 10 28 s _1 . Upper 
frame: Radial profile of band total Brightness vs. R (impact para- 
mater) for Phase angle = 180° for multiple Azimuthal angles. Pro- 
files of azimuthal angles (indicated by color coding in the online 
version) show no variation for this case (as should be expected) and 
overlap, appearing indistinguishable. Column density is included 
as the bold solid line (red in the online version) using a different 
y-axis scale. Lower frame: g-factors. Both the calculated g-factor 
(the higher group of lines) and g e ff, the observed brightness over 
column density, plotted with matching styles (colors in the online 
version) for each azimuthal angle (which also match those in the 
upper frame). 

breadth. This occurs primarily due to the temperature 
profile. We have used a fairly simplified profile, which 
can be scaled to a surface temperature parameter, but 
does not vary much otherwise between different model 
cases. This provides a straightforward "control" for this 
aspect of spectral change with aperture size. 

In the innermost coma near the nucleus, the temper- 
ature is quite warm (~100 - 200K), which leads to a 
broader band in the 10 km spectrum. The coma gas 
cools to a minimum (~20K) around 100-200 km out 
from the nucleus, which produces a much narrowed band. 
Since the Einstein A coefficient for the lowest J lines is 
higher than for the lines in the "wings" of the band, the 
cold temperature also increases the g-factor, even in opti- 
cally thick regions. At larger radii, the temperature rises 
again, but becomes less significant since the coma gets 
less dense and tends towards fluorescent equilibrium. Be- 
tween these regimes, in a "transition region," there are 
still optical depth effects, which can be more easily iso- 
lated as g-factors are less temperature controlled. 

Temperature is also a factor in determining doppler 



broadening and line width, which is proprtional to T 1 / 2 , 
so the ratio between line widths for the coldest and the 
warmest regions of the coma are about 2-3, for a given 
wavenumber. This may lead to temperature playing a 
significant role in the optical thickness of the coma to 
incident solar radiation. 

Temperature effects notwithstanding, the spectra from 
the denser near-nucleus regions of a coma show optical 
depth effects in several aspects. In addition to the total 
brightness no longer increasing linearly with production 
rate (and a corresponding reduction of g-factors), en- 
ergy is also dramatically shifted between lines within the 
band. 

The notable shifting of flux from R branch to P branch 
(evident in many of the thicker spectra), and to lower 
wavenumbers in both branches (as is most evident in the 
R ap =100 & 200 km spectra for Q = 10 28 s _1 ), are very 
noticeable optical depth effects. (See Sahai & Wannier 
1985, for an analytical discussion of similar effects.) 

This effect appears due to the branching ratio of a 
given pair of P and R branch lines originating in the same 
upper level, which generally (slightly) favors emission in 
the P branch line. In optically thick cases, repeated ab- 
sorption and emission of photons leads to a cumulative 
effect which favors the P branch over the R branch much 
more than in optically thin conditions (where it is prob- 
able that any emitted photon will not be re-absorbed 
before escaping the coma). 

Similarly, flux is "pushed" outwards in the branches, 
and more so in the P branch due to combination with the 
above effect. This is due to the lines closer to the center 
of the band becoming optically thick before those in the 
wings (both due to their higher Einstein coefficients and 
generally being more populated.) Flux initially emitted 
in lines that are optically thick will through repeated 
absorption and emission be forced out into lines that are 
less thick. 

3.5. The P/R Ratio: A Useful Heuristic of Optical 

Depth 

As seen above, the P branch total brightness and the R 
branch total brightness vary with respect to each other 
over different optical depths (as well as other factors, 
such as temperature distribution.) The ratio of the sums 
of P/R branches' brightnesses can be useful to alert an 
observer (or anyone analyzing observations) that they 
must in a given case beware of, and if possible account 
for, optical depth effects. 

This alone, would not be sufficient, as temperatures 
along a given line of sight are also a significant factor 
in controlling the P/R ratio, in the collisionally domi- 
nated inner coma. Colder population distributions will 
emit more in the lower lines (in both branches) for which 
the ratio of P /R for each pair of lines originating in the 
same upper state is greater. Also, dr will effectively vary 
inversely with linewidth, other factors being equal. 

Use of a model like ours can show where the P /R ratio 
is large due to temperature and where (its excess be- 
yond that value is) due to optical depth. In our optically 
thin, Qco = 10 26 s -1 , model the P/R ratio does not 
exceed ~1.44, even for aperture sizes dominated by the 
coldest portion of the coma. (Note, however, that this 
is an aperture averaged value. In Fig. 20 below, the 
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FlG. 17. — Aperture integrated spectra for Qco = 10 26 s 1 - Left side y-axis is aperture averaged brightness. Right y-axis is effective 
(line) g-factor (brightness/column density). Totals are indicated on each graph. 
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Fig. 18. — Aperture integrated spectra for Qco = 10 27 s 1 . Left side y-axis is aperture averaged brightness. Right y-axis is effective 
(line) g-factor (brightness/column density). Totals are indicated on each graph. 



peak value is slightly higher, ~1.5.) However, the ratio 
for corresponding aperture sizes in the Qco = 10 27 s _1 
and Qco = 10 28 s^ 1 cases is ~1.8 and ~2.4, respectively. 
Furthermore, in the thickest case modeled, even the spec- 



trum with aperture size of 2000 km has a ratio of ~1.4. 
Note that in all cases the 2 x 10 5 km aperture, which is 
dominated by the outer coma in fluorescent equilibrium, 
has a ratio of ~1.12. All of this indicates that a P/R 
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Fig. 19. — Aperture integrated spectra for Qco = 10 28 s 1 . Left side y-axis is aperture averaged brightness. Right y-axis is effective 
(line) g-factor (brightness/column density). Totals arc indicated on each graph. 



ratio in excess of ~1.4'- 
depth effects involved. 



1.5 is a warning sign of optical 



3.6. Further Discussion 



While it would be ideal to be able to derive a simple 
correction factor from the P/R ratio in such cases, alas, 
it is not exactly possible. However, a rough estimate of 
the degree of optical depth effects can be derived. 
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To do so, we create radial profiles of the P/R ratio, for 
both the observed emergent flux/brightness and the cal- 
culated value based on underlying populations wit h out 
attenuation of emergent light, as shown in Figs. [20l [2T1 
and 



CO 1-0 Band; Q= 1x1 ! V P/R Ratio vs. R 
(Phose angle = 90°, Multiple Azimuthal angles) 
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By cross-referencing the observed P /R ratiofor 
a given radial distance with the corresponding g-factor 
in Figs. [TUJ [TT] and [12] one can ascertain the "real" 
g-factor to use to calculate a correct column density from 
the observed flux. 



CO 1-0 Band; Q= 1x10"s"' P/R Ratio vs. R 
(Phase angle=90°, Multiple Azimuthal angles) 




Radius (km) 

Fig. 20.— For Q co = lQ 26 * -1 . Ratio of P branch vs. R 
branch total Brightness vs. R (impact paramater) for Phase angle 
= 90° and for multiple Azimuthal angles. Profiles of azimuthal an- 
gles show negligible variation for this case and overlap, appearing 
indistinguishable. (Profiles are indicated by color coding in the 
online version). 



CO 1-0 Band; Q= 1x10 ! V P/R Ratio vs. R 
(Phase angle = 90°, Multiple Azimuthal angles) 




Fig. 21.— For Qco = 10 27 s _1 . Ratio of P branch vs. R 
branch total Brightness vs. R (impact paramater) for Phase an- 
gle = 90° and for multiple Azimuthal angles. Profiles of azimuthal 
angles show minimal variation for this case and overlap, appearing 
nearly indistinguishable, except inwards of ^40-50 km. (Profiles 
are indicated by color coding in the online version). 



This heuristic is, of course, limited in use to to the 




Fig. 22.— For Q C o 



Ratio of P branch vs. 



branch total Brightness vs. R (impact paramater) for Phase an- 
gle = 90° and for multiple Azimuthal angles. Profiles of azimuthal 
angles show some variation in this case, but those for 0° , ±45° , 
and 90° overlap each other entirely in the calculated profiles and 
almost entirely in the observed. (Profiles are indicated by color 
coding in the online version). 

carbon monoxide X £ + band. Other spectra with P 
and R branches will have their own ratios, which can be 
derived by similar modeling. More complicated spectra 
may also, but such ratios would be more complicated to 
find than for cases with a simple two-branch structure. 

4. NEXT STEPS 

We plan to provide further useful results in forthcom- 
ing papers (currently in preparation) dealing with sim- 
ilar modeling to that presented here for H 2 and CO2, 
as well as a more in depth treatment of CO alongside 
those species. We will present model results including 
morphology and comparison with the in-situ spectral ob- 
servations of the Deep Impact and EPOXI missions. The 
model and code have already been implemented to pro- 
duce those results. 

Further work, not yet implemented, will include sev- 
eral planned improvements to our model. A more accu- 
rate treatment of radial velocites and doppler shifts in 
lines due to them is a highly desirable improvement. (At 
present, only a thermal doppler profile is used for cal- 
culating absorption of solar radiation.) We neglected a 
precise treatment until now primarily due to the nature 
of our originally intended problem. The near nucleus 
morphology which we intended to model is represented 
by cones expanding radially in one direction, a geometry 
that is likely to reduce doppler effects. However, when 
modeling the whole coma as a sphere, as in the results 
presented here, this assumption is no longer valid. Fur- 
thermore, in as much as we are primarily looking at the 
overall band shape as opposed to individual line shapes, 
the effect of doppler shifts on the spectra is expected to 
be less than if one were modeling specific lines (as is often 
the case for other spectral regimes). However, the effect 
of this neglect would be to increase the optical thickness, 
and therefore our results may be over-estimating optical 
depth effects. If so, the effects we describe would still be 
observable, but for higher production rates than those 
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modeled. 

In addition to more accurate radial velocities, more 
flexible radial temperature and density profiles are 
planned, so as to be able to model deviations from a 
very simple Haser model, (e.g. Volatiles produced from 
icy grains or large chunks and not solely from the nu- 
cleus' surface, as per A'Hearn, et al. 2011, and/or pho- 
todissociation of molecules, and corresponding creation 
of daughter species.) 

Computational limits are currently a limiting factor in 
how optically thick and how refined (in terms of granu- 
larity of conditions, in that more variation requires more 
regions) the modeled cases can be. As of now, the maxi- 
mum production rates we can deal with are on the order 
of Q — 10 28 s _1 , and somewhat less for CO2 than the 
other molecules. We are planning to address these con- 
cerns with algorithmic improvements. Running the code 
on faster and more powerful computers is also a possibil- 
ity. 

5. CONCLUSIONS 

We have demonstrated our model's usefulness in un- 
derstanding emission spectra of cometary comae. There 
are several possible effects of optical depth that could 
lead observers to mistaken conclusions regarding the cal- 
culated abundances, or other characteristics, of species of 
interest. The moral of the story: Ignore radiative trans- 
fer and optical depth effects at your own peril! 

Although designed specifically with comets in mind, 
our model and code are versatile enough to be used in 
other radiative transfer problems as well. Parameters 
that define a specific comet model or other problem, in- 
cluding molecule of interest, size of nucleus and radial 
shells, production rate, morphology (if any), incident ra- 
diation, etc. are all fairly flexible. Thus our adaptation 
of Coupled Escape Probability to an asymmetric spheri- 
cal situation has created a very useful tool for modeling 
cometary spectra, as well as other spherical astrophysical 
phenomena. 



APPENDIX 

ALGORITHM IMPLEMENTATION & TECHNICAL 
DETAILS 

We have coded the aglorithm described above in the 
CH — h language, using numerous functions from Press, et 
al. 1992, primarily to implement numerical integration of 
functions (with odeint and stifbs and associated func- 
tions) and solution of N-dimensional non-linear matri- 
ces with Newton's Method (using newt, and associated 
functions). The bulk of the coding, which implements 
the radiative transfer algorithm in spherical geometry is 
our own. 

A major practical limitation of our algorithm is the 
matrix size; since Newton's method requires (repeated) 
0(N 3 ) matrix solving operations (for which newt uses 
the brute force approach of ludemp and lubksb), the al- 
gorithm can get prohibitively slow for large matrices. For 
example, on an Intel Core computer (with CPU speed 
of 2.9 GHz and 7.6 GB of memory) running Scientific 
LINUX 6.3 (Carbon), when N > 15,000, a solution may 



take one or more days. For greater sizes, even a week. 
The matrix size equals the number of molecular levels 
used times the number of regions. The molecule (and 
band/s) being modeled determines the first value, with 
the latter value demanding to be increased with greater 
optical depths and production rates. Depending on the 
species of interest, the maximum practical production 
rates we can currently manage on such a system are of 
O(10 2S )s~ 1 . The most simple workaround is to use a 
more powerful computer. Algorithms for solving sparse 
matrices faster than the above functions can also be used, 
and we have begun to explore this option. 

In implementing the CEP method, and our adapta- 
tion to spherical geometry, we have created C++ classes 
for diatomic and triatomic molecules. The object ori- 
ented programming style of C++ lends itself ideally to 
being able to switch molecules easily. The Diatomic 
molecule class and its subclasses (currently implemented 
for CO and SiO) calculate energies and Einstein coef- 
ficients based on constants (taken primarily from Kru- 
penie 1966, for CO) that are included in the code. The 
Triatomic class, which can actually be used for other 
polyatomic molecules as well (or the aforementioned di- 
atomic molecules themselves, for that matter), must be 
provided with energies and coefficients from some out- 
side source in formatted input files. We used the HI- 
TRAN database (see Rothman, et al. 1998) to supply 
these values for CO2 and H 2 0. This approach has the 
versatility to handle many other molecules with a mini- 
mal effort of "data massaging" to get the data into the 
proper format. We have also created a large number of 
classes which encapsulate the spherical geometry and the 
attendant calculations. These are all "controlled" by the 
Comet class that reads parameters for a given case from 
a "comet definition file" (a text file of key/value pairs, 
essentially) and, using the above classes, sets up and runs 
the model for a given case. 

Although designed specifically with comets in mind, 
our code is versatile enough to be used in other spheri- 
cal radiative transfer problems as well. Straightforward 
input files describe all the required and optional parame- 
ters for a specific comet model or other problem, includ- 
ing molecule of interest, size of nucleus and radial shells, 
production rate, morphology (if any), incident radiation, 
etc. 

The C++ code outputs a file containing a point-by- 
point line-by-line spectral mapping, as described above, 
for one or more specified viewing orientations. This data 
can then be presented in multiple formats as described 
above in Section 3. We have implemented this data pre- 
sentation portion of the model using various short IDL 
programs which we have developed, each specifically for 
the purpose of producing a given type of graphical out- 
put. 
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